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We use a long, all-atom molecular dynamics (MD) simulation combined with theoretical modeling 
to investigate the dynamics of selected lipid atoms and lipid molecules in a hydrated diyristoyl- 
phosphatidylcholine (DMPC) lipid bilayer. From the analysis of a 0.1 jjls MD trajectory we find 
that the time evolution of the mean square displacement, ([£r(£)] 2 ), of lipid atoms and molecules 
exhibits three well separated dynamical regions: (i) ballistic, with {[Sr(t)} 2 ) ~ t 2 for t < 10 fs; (ii) 
subdiffusive, with ([Sr(t)] 2 ) ~ t^ with /3 < 1, for 10 ps < t < 10 ns; and (hi) Fickian diffusion, with 
{[Sr(t)} 2 ) ~ t for t > 30 ns. We propose a memory function approach for calculating {[Sr(t)] 2 ) over 
the entire time range extending from the ballistic to the Fickian diffusion regimes. The results are 
in very good agreement with the ones from the MD simulations. We also examine the implications 
of the presence of the subdiffusive dynamics of lipids on the self-intermediate scattering function 
and the incoherent dynamics structure factor measured in neutron scattering experiments. 
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I. INTRODUCTION 

Lipids are polymer molecules composed of hydrophobic 
acyl chains attached to a hydrophilic polar head group. 
In the presence of a polar solvent (e.g., water) lipids can 
spontaneously self-assemble to form a bilayer membrane 
(Fig. [T]). The fluid (L a ) phase of lipid bilayers behaves 
as a two dimensional (2D) fluid and the lateral 2D self- 
diffusion coefficient of individual lipids within a leaflet 
of the bilayer has been determined both experimentally 
by a variety of methods [IJI2l3[l[3[l[3ia[3[l0] 
and through computer simulations [TTJ [12j [13] [14j [15] 
EE El [THJ [19] . The experiments suggests that there are 
at least two relevant length/time scales associated with 
the lateral diffusion of lipids in a bilayer. Experiments 
designed to probe motion on picosecond (ps) time scales 
measure a diffusion coefficient D\ that can be one to two 
orders of magnitude larger than the diffusion coefficient 
D2 measured in experiments which probe motion on tens 
or hundreds of nanosecond (ns) time scales. 

There have been several models used to explain the 
difference in the diffusion coefficients D\ and D2. Vaz 
and Almeida [20] suggested a free- volume jump-diffusion 
model where the lipid moves in discrete steps when a 
void forms next to the lipid. After a jump, the lipid can 
either return to its original position or another lipid can 
jump into the empty space left behind. In this model, the 
short time diffusion coefficient D\ is associated with the 
"rattling" motion of the lipid inside the "cage" created 
by the surrounding lipids. The jump diffusion mecha- 
nism was investigated by Falck et al. [TT] in a simulation 
of diyristoyl-phosphatidylcholine (DMPC) bilayers, and 
they did not observe enough jump events to provide ev- 
idence for the jump diffusion model. To explain their 
MD simulation results of a DMPC bilayer, Wohlert and 
Edholm [12] proposed a model where the motion of the 
lipid in the plane of the bilayer can be regarded as diffu- 
sion (with a coefficient D\) confined in a circular "cage" 



whose center undergoes free diffusion (with a coefficient 
D2). The two diffusion coefficients obtained by fitting 
their simulated data to the analytical solution of this 
model were in reasonable agreement with the experimen- 
tal values. 

None of the above diffusion models explicitly take into 
account the microscopic polymeric structure of lipids, but 
effects of this structure has been observed in simulations. 
For example, Doxastakis et al. [13 observed in a MD sim- 
ulation of l,2-Dipalmitoyl-sn-Glycero-3-Phosphocholine 
(DPPC) bilayers that the atoms at the ends of the lipid 
tails fluctuate more than those in the head group, but the 
overall diffusion of the lipid is limited by the diffusion of 
the head group. This observation led them to examine a 
diffusion in a sphere model where the size of the sphere 
depends on the position of the chain. However, the au- 
thors describe this model as "not particularly accurate" , 
and the model was extended to include a distribution of 
spheres at each tail position. 

The dynamics of the atoms in lipid molecules are more 
complex than those in ordinary liquids. In a simple liq- 
uid [21] atoms move ballistically at short times. Thus the 
mean square displacement (MSD) {[Sr(t)] 2 ) « ((vt) 2 ) ~ 
t 2 , which is followed by a crossover to Fickian diffusion, 
characterized by ([Sr(t)] 2 ) ~ t for long times. In dense 
fluids a caging effect, where the atoms are trapped by 
their neighbors, is observed between the ballistic and dif- 
fusive regimes, leading to a plateau in ([Sr(t)] 2 ). In a lipid 
bilayer, the motion of lipid atoms is further complicated 
by the polymeric structure, characterized by connectivity 
and flexibility, of the lipids. 

A more realistic description of the diffusion of lipids 
should be based on theoretical models designed for poly- 
mers. Two such theories are the Rouse model [22] and 
the mode-coupling theory [23] [24] . In the Rouse model 
over-damped Brownian motion is assumed for the in- 
dividual monomers with a harmonic potential connect- 
ing adjacent monomers. The motion is subdiffusive with 
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([Sr(t)} 2 ) ~ t 1 ! 2 then there is a crossover to diffusive mo- 
tion ([Sr(t)] 2 ) ~ t 1 at later times. The mode-coupling 
theory (MCT), generally used to study the dynamics of 
glass- forming liquids, has recently been extended for flex- 
ible macromolecules and polymers [23j [24] , though it has 
not been applied to study lipid bilayers to our knowledge. 
The MCT is a method used to approximately calculate 
space-time correlation functions (e.g., intermediate scat- 
tering function) by solving a set of integro-differential 
equations obtained from first principles. An important 
feature of the extended polymer version of the mode 
coupling theory is the prediction of a sub-diffusive re- 
gion between the short time ballistic and the long time 
Fickian diffusion regimes. Specifically, the theory pre- 
dicts ([Sr(t)] 2 ) ~ , with (3 typically between 0.5 and 
0.65. Besides other parameters, the exponent f3 de- 
pends on the length of the polymer chain and approaches 
0.5, the Rouse limit, for an infinite chain of identical 
monomers [23, 24 . Since the only input to the theory is 
the static structure, the subdiffusion predicted by the ex- 
tended MCT is due to the polymeric structure and not to 
the specifics of the interactions. Thus, one expects that 
the dynamics of the atoms in a lipid bilayer should ex- 
hibit a pronounced sub-diffusive regime before it crosses 
over to normal Fickian diffusion. We argue that the ex- 
tended sub-diffusive region in the MSD is the cause of the 
difference between the experimentally measured diffusion 
coefficients D\ and D^. 

In this paper, we report an extensive computational 
investigation of the complex dynamics of both individual 
lipid atoms and entire lipid molecules in the fluid phase 
of a hydrated DMPC bilayer by employing a microscopic 
description of the system. Our study is based on a 0.1 /is 
long all atom molecular dynamics (MD) simulation of 
the DMPC bilayer. By calculating the time evolution of 
the mean square displacement of lipid atoms in the plane 
of the bilayer for a time interval ranging from 10 _15 s to 
10 _7 s, we identify the short time (t < 0.1 ps) ballistic re- 
gion, the long time (t > 30 ns) Fickian diffusion region, 
and in between an extended subdiffusive region character- 
ized by a power law {[Sr(t)] 2 ) ~ t@, with (3 < 1. We find 
that (3 depends on the atom type and, most importantly, 
on its position within the lipid molecule. The dynam- 
ics is different for atoms in the lipid tail versus atoms 
in the head group. We also examine the implications of 
the heterogeneous subdiffusion of the individual atoms 
on the diffusion of the lipid molecule as a whole, and on 
the interpretation of neutron scattering experiments. 

This paper is organized as follows. In Sec. [TT] we de- 
scribe our MD study of a solvated DMPC lipid bilayer. 



In Sec.yTywe provide a detailed numerical and theoreti- 
cal study of the mean square displacement of lipid atoms 



and lipid molecules. In Sec. IV we examine the implica- 
tions of this study for the incoherent intermediate scat- 
tering function, and the dynamic structure factor that is 
measured in neutron scattering experiments. Lastly, in 
Sec. [V] we present a discussion and our conclusions. 
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FIG. 1: (color online) (a) The simulated DMPC lipid bilayer 
system, (b) A DMPC lipid molecule. The atoms whose dy- 
namics is investigated in this work (P, CH and CT) are high- 
lighted as black spheres. Figures created using VMD [25] . 



II. MD SIMULATION OF DMPC BILAYER 

To investigate the dynamics of selected lipid atoms 
and of entire lipid molecules, we performed a 0.1 fis long 
all atom MD simulation of a fully solvated DMPC lipid 
bilayer (Fig. [I]). The pre-equilibrated structure of the 
bilayer was obtained from Mikko Karttunen's web site 
(www.apmaths.uwo.ca/ ^mkarttu/downloads.shtml) [26] 
and contained 128 lipid molecules. The system was sol- 
vated by adding two 11 A thick layers of water to each 
side of the membrane using the Solvate plugin in VMD 
[25] . The final system contained a total of 2577 TIP3 [27] 
water molecules. The MD simulation was performed with 
NAMD-2.6 [28] using the CHARMM27 [29] force field. 
The equations of motion were integrated with a multiple 
time step algorithm with time steps of 1 fs for bonding 
interaction, 2 fs for non-bonding interactions, and 4 fs for 
long-range electrostatic interactions. The non-bonded in- 
teractions were cutoff at 12 A with a smooth switching 
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function starting at 10 A. Long-range electrostatic inter- 
actions were computed by employing the smooth Particle 
Mesh Ewald method [30] with a grid spacing of 1 A. 

The system was brought to T = 303 K and normal 
pressure (p = 1 atm) through several stages of equilibra- 
tion. First, the system was subjected to 6 x 10 4 energy 
minimization steps by harmonically restraining the phos- 
phorous (P) atoms of the lipid headgroups along the nor- 
mal direction to the surface of the membrane. Next, the 
system was gradually heated to the desired temperature 
T = 303 K during a 0.5 ns period. After removing all the 
restraints, the system was equilibrated through a 75 ns 
long NpT simulation. At the end of this process the area 
per lipid (APL) was 56.2 A 2 . To increase the APL to 
the desired value of ~ 60 A 2 , we performed NVT sim- 
ulations and gradually increased the size of the system 
in the x?/-plane of the membrane. The final system size 
was 62 x 62 x 58.5 A 3 , with an APL of 60.06 A 2 . The 
equilibration process was concluded with an additional 
10 ns NVT simulation. 

A production run of 0.1 /is was performed in the NVT 
ensemble. We employed a Langevin thermostat with a 
coupling constant of 0.05 ps _1 . The coordinates of all 
the atoms were stored every 2 fs for the first 100 ps, then 
saved every 20 fs for the next 10 ns, then every 100 fs 
for the next 90 ns. This allowed us to study the short, 
intermediate, and long time dynamics of the individual 
atoms and lipids. A snapshot of the lipid bilayer system is 
shown in Fig. [lji. The MD simulations were carried out 
on 40 CPUs of a dual core 2.8GHz Intel Xeon EM64T 
cluster with a performance of around 0.2 days/ns. 



III. MEAN SQUARE DISPLACEMENT 

In this section we examine the time dependence of 
the lateral mean square displacement (MSD) of selected 
atoms within a lipid, and of entire lipids within the bi- 
layer. We compare the results with the general predic- 
tions of the extended MCT for polymers [23j El] , and we 
propose a memory function method capable of describing 
the time dependence of MSD in lipid bilayers. 

We begin by outlining some of the qualitative predic- 
tions of the extended MCT for polymers. To derive the 
extended mode coupling equations, one first applies the 
Mori-Zwanzig projection operator formalism in the so 
called site-site formalism. This leads to a set of cou- 
pled integro-differential equations for monomer density 
fluctuations. The MCT is a set of approximations of the 
memory kernel that permits the corresponding integro- 
differential equations to be numerically solved with the 
only input being the static structure. The monomer and 
the polymer center of mass MSD can be computed from 
the solution of the mode-coupling equations from the 
small wave- vector, limit (see Ref. l24l for details). For 
the MSD ([Sr a (t)} 2 ) = ([r a (t) - f a (0)] 2 ) for monomer 
a, the MCT predicts an initial ballistic region where 
([5r a (t)] 2 } ~ t 2 . After the ballistic regime, the growth 
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FIG. 2: Lateral mean square displacement of the P (squares), 
CH (triangles) and CT (open circles) atoms, along with the 
center of mass (CM) of the lipids (closed circles). The se- 
lected atoms (P, CH and CT) are highlighted in Fig.[l] Inset: 
a(t) = d\n([Sr CM (t)] 2 )/d\n(t). The dashed horizontal lines 
correspond to the ballistic (a = 2), sub-diffusive (a = 0.677) 
and Fickian diffusion (a = 1) values. 



of the mean square displacement is suppressed due to a 
local caging effect caused by the trapping of the polymer 
by its neighbors. At larger times there is a crossover to 
Fickian diffusion with ([Sr a (t)} 2 } = 2dD L t, where d is 
the dimension of the system and Dl is the long time self 
diffusion coefficient for the monomer. In this crossover re- 
gion there is a sub-diffusive region where ([5r a (t)] 2 ) ~ t@ 
with (3 < 1. The value of (3 is smaller for longer chains 
and has a limiting value of 0.5 for an infinitely long chain. 
Note that the sub-diffusive region is absent for simple liq- 
uids, and is a polymer specific feature. 

To test the MCT predictions for the lipids, we exam- 
ined the motion of three representative atoms in the lipid 
(see Fig. [l]3), namely the phosphorous (P), the carbon 
C15 (CH) atoms in the lipid head group, and the carbon 
atom C214 (CT) located at the end of one of the lipid 
tails. Because the surroundings of the head group atoms 
is different from those in the lipid tails, we expect that 
the MSD of the P and CH atoms will show distinctive fea- 
tures from that of the CT atoms. Shown in Fig. [2] is the 
time evolution of the lateral MSD of the P atom (squares) 
([Sr p (t)] 2 ) = {\r p (t) — r p (0)| 2 ) on a log- log scale, which 
follows the qualitative predictions of the MCT. Specifi- 
cally, an initial ballistic region is followed by a narrow 
caging region that gradually crossover to linear diffusion 
through an extended sub-diffusive region characterized 
by ([5r p (t)} 2 ) ~tP, with [3 <1. 

Also plotted in Fig. [2] is the MSD for the CH (triangle) 
and CT (circles) atoms. The ballistic, the sub-diffusive, 
and the diffusive time scales are identifiable for all the 
atoms. The MSD in the ballistic regime is the same for 
CH and CT since they have the same mass and are at the 
same temperature, but their mean square displacements 



4 



differ starting in the caging region. The MSD for CT is 
larger than that of CH, which indicates that the size of 
the cage is larger for atoms at the end of the tail than 
atoms close to the head group. This difference is a direct 
consequence of the structure of the lipid bilayer, and has 
been observed in previous simulations [I2j [13] . 

For all the atoms the crossover to Fickian diffusion 
begins around 10 ns and ([8r a (t)] 2 ) « 100 A 2 , which ap- 
proximately corresponds to the nearest neighbor distance 
of the lipids 10 A). Thus the linear diffusion can be 
observed only after the lipid has moved around one lipid 
diameter. 

To determine the power-law exponent (3 a for atom a 
in the sub-diffusive region, we calculated 

n ,a (f , dM([5r*m} t d([Sr%t)} 2 ) 

[ } d\n(t) ([Sr a (t)} 2 ) dt ' 1 ; 

and averaged the result from t = 0.01 ns to t = 10 ns, 
the region where a a (t) is approximately constant. The 
result for a a (t) for the center of mass (CM) of the lipids 
is shown in the inset to Fig. [2] The horizontal lines in 
the inset correspond to the ballistic {f3 = 2), sub-diffusive 
(/3 = 0.677) and Fickian diffusion (/3 = 1) regimes. The 
value of f3 a depends on the atom type and ranged from 
0.677 ± 0.005 for the CM to 0.427 ± 0.004 for CT, the 
carbon atom at the end of the tail. The dynamics in 
the ps time scale strongly depends on the position of the 
atom within the lipid. The atoms at the end of the lipid 
tails have more freedom to move, and hence a larger cage. 
However, these atoms are connected to those within the 
head group, which ultimately determine their diffusion 
on long time scales. 

For t > 30 ns, a a (t) « 1 for each of the individ- 
ual atoms (a = P, CH, CT) and the center of mass, 
thus the lipids appear to be undergoing Fickian dif- 
fusion. We calculated the diffusion coefficient D a L = 
lim £ ^ 00 ([(5r a (t)] 2 )/4t by fitting ([Sr a (t)} 2 ) for t > 30 ns 
to a straight line. The diffusion coefficients (Dz,-msd) are 
shown in Table ffl and their values range between D2 T = 
1.67 x 10~ 7 cniV 1 and L>£ H = 1.27 x 10~ 7 cmV 1 . 
These values are within statistical error and are close to 
the diffusion coefficient of the CM. However, these diffu- 
sion coefficients are almost a factor of two larger than the 
value of L> exp = 0.69 x 10~ 7 cnrV 1 obtained from FRAP 
measurement [10] of DMPC at 34° C. This difference may 
be due to a combination of several factors, e.g., imperfec- 
tion of the empirical CHARMM27 force field used in the 
MD simulations, discrepancy between the APL of the ex- 
perimental and the simulated systems, finite size of the 
DMPC bilayer system and insufficient sampling. How- 
ever, the fact that we have obtained similar long time 
diffusion coefficient for all lipid atoms indicates that our 
MD trajectory is sufficiently long to capture the linear 
diffusion regime. 

We propose an approach based on the Zwanzig- 
Mori projection operator method to model the MSD, 
([5r a (t)] 2 ), for a lipid atom a in terms of five fitting 
parameters. We start from the equation of motion for 



TABLE I: The diffusion coefficient calculated from a linear fit 
to the the mean square displacement and the memory function 
approach described by Eqs. [3] and [4] for the atoms shown in 
Fig. [l] and the center of mass (CM). The experimental value 
obtain from FRAP experiments is 0.69xl0 -7 cm 2 s _1 [TO] . 
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the density, p a (q : t), autocorrelation function ^(q : t) = 
(Pa(—Q,Q)Pa(q,t)) of a tagged atom a at wave-vector q 

m 

dUsil, t)+q 2 v 2 a r s (q, t)+ f M a (q, t-s)d s <j> a s (q, s)ds = 0. 

Jo 

(2) 

Here v a = y^ksT '/ra a , ks is the Boltzmann's constant, 
T is the temperature, m a is the mass of atom a, and 
M a (g, t) is the memory kernel. In general, MCT refers to 
an approximation for M a (q,t). The equation of motion 
for the MSD in 2D can be obtained from ([5r a (t)] 2 ) = 
41im^ [l - ^ste^l/V which gives 

dt([Sr a (t)} 2 ) + f M$(t - s)([Sr a (s)¥)ds = 4v% (3) 

where Mg(t) = Hm g _ M a (q, t). 

We propose the following ansatz for the memory func- 
tion 

M S(t) = s { t)M + J fg^, (4) 

where S(t) is the Dirac delta function. The physical 
meaning of the fitting parameters r§ , B ai rf , and 
f3 a are clarified next. 

If B = 0, then Eqs. describe the free diffusion 

of a Brownian particle with a diffusion coefficient D3 = 
k B T/(r 3 m). Indeed, when 5 = 0, Eq. ^ can be solved 
analytically and ([Sr(t)} 2 ) = 4D 3 t - 4D 3 r 3 (l - e' 1 '^). 
For t <C r 3 , ([Sr(t)} 2 } w (2ksT /m)t 2 corresponding to 
the ballistic regime, while for t ^> r 3 ([Sr(t)] 2 ) w 4D 3 1 
corresponding to the Fickian diffusion regime. To illus- 
trate the crossover between these two limiting cases we 
plotted a(t) defined by Eq. 0, Fig. [3] For free diffusion 
a(t) changes smoothly from a = 2 (ballistic dynamics) to 
a = 1 (Fickian diffusion) as shown in Fig. [3] (dot-dashed 
line) for r 3 = 10 fs. 

By neglecting the power-law term in the memory 
kernel and setting M°(t) = S(t)/r 3 + Bexp(-t/ri), 
([Sr(t)] 2 ) can be calculated analytically, though the ex- 
pression is too cumbersome to be shown here. In this 
case, first the MSD crosses over from the ballistic regime 
(for t <C T3) to a caging region (for r 3 <C t <C T\) charac- 
terized by a plateau in ([Sr(t)] 2 ) (with a = 0) as shown 
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TABLE II: Fit parameters of the mean square displacement 
using the memory function approach, Eqs. ([3} and Q. 
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FIG. 3: a(t) for MSDs calculated from Eq. Q for different 
memory kernels M(t) as described in the text. 



in Fig. [3] (dashed curve). Then there is a crossover from 
the caging to the Fickian regimes that takes place around 
r\. The extent of the caging region is determined by the 
difference between the two time scales T3 and t\. 

The power-law term in the memory kernel, Eq. Q, is 
responsible for the anomalous subdiffusion in the MSD. 
To illustrate we have numerically calculated the MSD 
and the corresponding a(t) for M°(t) = S(t)/r 3 + B/[1 + 
(t/ r 2)~^]- As shown in Fig. [3] (long-dashed curve), in 
a broad time interval centered around T2 = 0.1 ps and 
extending over four decades, the dynamics slowly crosses 
over from ballistic to sub-diffusive (with (3 < 1) behavior. 

To capture all the features of the mean square displace- 
ment of lipid atoms in a bilayer membrane one needs to 
consider the full kernel Q with all three time scales. A 
representative result for a(t) obtained with the proposed 
memory function method is shown in Fig. [3] (solid line). 

Based on the above analysis, one can identify n, T2, 
and T3 as the characteristic time scales for the crossover 
from the subdiffusion to Fickian diffusion region, the on- 
set of the subdiffusion regime, and the crossover from 
ballistic to caging region, respectively. To determine the 
times Ti, T2, and T3 we performed a least squares fit us- 
ing Eqs. ([3} and @ to ([Sr a (t)} 2 ). The fits are shown 
in Fig. [4] and the fit parameters are given in Table [TTJ 
The inset to the figures show the corresponding a a (t) 
for the fits and calculated from the MD simulation. For 
all the individual atoms there are features in a a (t) for 
t < 1 ps that are not captured by our memory function 
method. The crossover to Fickian diffusion also appears 
to be sharper in the simulations than predicted by the 
theory. However, the memory function theory fits the 
MSD from MD simulation well for over eight decades in 
time. 

The three time scales are separated by at least an or- 
der of magnitude, and define three clearly separated dy- 
namical regimes of the lipid bilayer system. The short 
time scale T3 ~ 10 fs characterizes the crossover between 
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the ballistic and caging regions. The intermediate time 
scale T2 ~ 1 ps defines the onset of the sub-diffusive re- 
gion. Finally, the long time scale t\ ~ 10 ns identifies 
the crossover from subdiffusion to Fickian diffusion. 

We conclude this section by noting that the memory 
function approach provides another means to calculate 
the diffusion coefficient Dl. Inserting ([Sr(t)] 2 ) « 4Dt 
(valid in the t — > 00 Fickian diffusion regime) into Eq. ([3| , 
one finds that 



1 -1 



Mg(t)dt 



(5) 



The diffusion coefficient Dl should be independent of the 
atom used in the calculations. The diffusion coefficients 
(.D^-theory) for the P, CH, CT atoms and the center of 
mass calculated using Eq. ([5| are given in Table [TJ These 
values are in good agreement with those determined di- 
rectly from the MD simulation, but are systematically 
smaller. Furthermore, like the fits to the asymptotic be- 
havior of the MSD, they are nearly equal. Note that the 
fitting parameters differ by an order of magnitude for 
different atoms in the lipid, but the long time diffusion 
coefficient is the same for all the atoms. 



IV. INCOHERENT INTERMEDIATE 
SCATTERING FUNCTIONS 



To understand the implications of the large sub- 
diffusive region in the analysis of neutron scattering ex- 
periments, we examined the self (incoherent) intermedi- 
ate scattering function (SISF) for different atoms in the 
lipid and the lipid center of mass. We start this section 
with some theoretical background that relates the mo- 
ments of atomic displacements with the SISF in terms 
of a cumulant expansion. Next, we examine the viabil- 
ity of the second (Gaussian) and fourth order cumulant 
approximations for calculating the SISF for the selected 



atoms studied in Sec. Ill This will allow us to quantita- 
tively relate the MSD to the SISF and the interpretation 
of neutron scattering data. Finally, we perform the same 
analysis for the lipid hydrogens, which dominate the in- 
coherent neutron signal. 
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t [ ns ] t [ ns ] 

FIG. 4: The fits to the mean square displacement using the memory function approach, Eqs.[3]and^] The solid line is (Sr 2 (t)) 
calculated from the simulations and the dashed lines are the fits. The insets show the exponent a(t) — <91n([(!)r a (£)] 2 )/<91n(£) 
as a function of time. A running average was used to decrease the noise in the simulation results. 



A. Background and Theory 

The differential cross section of the quasi-elastic scat- 
tering of neutrons in a solid angle dQ with an energy 
transfer Tiuj can be expressed as [31] 



oc 



J2 (&SJ 2 S?(q,u) + Y, h coh h co h ST h {q,u), 



da 2 
dQdco 

n n,m 

(6) 

where n, m are atom-type indices, while bf nc [b™ oh ] and 
S™{q,uS) [S™%(q, uj)] are, respectively, the incoherent [co- 
herent] scattering length and dynamic structure factor. 
S™(q,uj) contains information about the single particle 
motion of nuclei of type n, and in principle it can be 
used to determine the self diffusion coefficient Dl- 

In a computer simulation, S™(q,u) is obtained from 
the Fourier transform of the SISF 



\m=l 



;(o)- 
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(7) 



where the summation goes over all N n atoms of type n 
in the system. Note that the only difference between 
I™(q,t) and <^ (#,£), discussed in Sec. |III| is that I™(q,t) 
is defined for a subset of atoms with the same scattering 
length, while 4^(q,t) is defined for a single tagged atom. 
The fourth order cumulant expansion of Eq. ^ reads 



:([8r n (t)] 2 ) 



1 , i ( ?([*"<!)?) ) ^ 



7 2 "(*) 



(8) 



where 



7?(*) 



2{[Sr n (t)] 2 ) 



1. 



(9) 



Equations ^ and ([9| express the SISF is terms of the 
MSD and the non-Gaussian parameter (t). If the 
distribution of lateral displacements of the lipid atoms 
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G"(r, t) = (<5(r - \ f n (t) - r"(0)|)) is Gaussian in space, 
then 



I£(q,t) =exp -4 



o([<5r"(i)] 2 ) 



(10) 



For 2D Fickian diffusion ([Sr n {t)} 2 ) = 4D L t, G"{r,t) is 
Gaussian in space, and 



I?(q,t)=e 



_ c -q 2 D L t 



(11) 



Note that, while the distribution function of the displace- 
ments of atoms that undergo Fickian diffusion is Gaus- 
sian (having a width that increases as the square root of 
time), the converse in not true in general. 

The dynamic structure factor corresponding to the 
SISF given by Eq. (11) is a Lorentzian [21] 



q 2 D L 



1 



7T (q 2 D L ) 2 



(12) 



which is the basis of evaluating the diffusion coefficient 
Dl from neutron scattering experiments. To determine 
the diffusion coefficient, it is generally assumed that the 
translational motion of the lipid is decoupled from other 
motions of the atoms within the lipid (e.g., vibrations, 
rotations, conformational changes, etc.). The resulting 
S™(q,cu) is a convolution of these other motions with 
Eq. ( fl2| ). From the above discussion, however, it should 
be clear that the reliable experimental determination of 
the diffusion coefficient Dl based on Eq. (12) is possible 



only when the MSD of the lipid atoms and/ or the lipid 
center of mass increases linearly in time. 



B. Individual Atoms 

In this sub-section we calculate the SISF for the P, 
CH and CT atoms in the lipid, and examine the validity 
of the expansion given by Eq. (|8| for these atoms. This 
analysis will give insight into the different dynamics along 
the lipid and aid in the analysis of the dynamics of the 
hydrogen atoms discussed in the next sub-section. 

Shown as symbols in Fig. [5] is I™(q,t) for q G 
{1.42,0.75,0.5} A" 1 calculated for the selected P, CH 
and CT atoms. The peak in the static structure factor for 
the lipid tails (head groups) occurs around q — 1.42 A -1 
(q = 0.75A -1 ). The Gaussian approximation (solid 
lines), Eq. ( [To] ) , of the SISF appears to be excellent for 
the P atoms(open squares) for all q values considered, 
but noticeable deviations are present for the carbon (CH 
and CT) atoms. Note, however, that as q decreases the 
Gaussian approximation becomes increasingly better for 
the carbon atoms as well. 

In order to quantify how well the distribution function 
(r, t) of the lateral displacements of lipid atoms of type 
n can be approximated by a Gaussian, we calculate the 
non-Gaussian parameter (t) given in Eq. (|9|. Shown 
in Fig. [6] are j%(t) for the atoms P, CH, CT, for all the 
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FIG. 5: Self intermediate scattering function for P (squares), 
CH (triangles) and CT (open circles) atoms for q = 1.42A -1 
(top), q = 0.75A -1 (middle) and q = 0.5A -1 (bottom). 
The solid curves correspond to the Gaussian approximation 
Eq. (PLOT). 



carbon atoms within the lipid, and for the center of mass 
of the lipids. For the carbon atoms, there is a peak in 
72 (t) around the beginning of the sub-diffusive region in 
the MSD, (i.e., t « 10 ps). The peak heights for CH 
and CT atoms are at least a factor of three smaller than 
when the displacements are averaged over all the carbon 
atoms. For times when 72 (t) is small, Eq. (10) is a good 
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FIG. 6: Time dependence of the non-Gaussian parameter 
72 (£) f° r the P (square), CT (circle) and CH (triangle) atoms, 
as well as, averaged over all the carbon atoms within the lipid 
(diamond), and calculated for the center of mass (solid line). 



approximation, which can be seen by comparing Figs. [5] 
and[6j Furthermore, it is evident that the approximation 
given by Eq. ( 10 ) is very poor when applied to the average 



of the displacements of all the carbon atoms within the 
lipid. Notice that 72 (t) ~ when t « 30 ns for all n, 
which is the time scale for the crossover to the linear 



region in the MSD (t\ in Table |II|) found in Sec. Ill 
We examined the approximation given by Eq. 



CT, the atom with the largest value of 72 (t). Note that 
for short and long times G™(r,t) is approximately Gaus- 
sian is space, which is evident by the small values of 
the non-Gaussian parameter 7 2 ~ 7T (£) for these times (see 
Fig. [6|. For intermediate times that correspond to the 
subdiffusion dynamics, 72 7T (^) is finite (especially for the 
carbon atoms) and, according to Eq. (J8|, one expects 
noticeable deviation from the Gaussian approximation, 
Eq. (10). A comparison between the direct calculation, 
the Gaussian approximation given by Eq. jT0| ), and the 
first non- Gaussian correction given by Eq. ^ of the SISF 
for the CT atoms is shown in Fig. [7] for two scattering 
vectors qi = 1.42 A -1 and q 2 = 0.75 A -1 . The figure 
shows that the first correction term matches very well 
the direct calculation of I™(q,t) for all t and both q val- 
ues, while the Gaussian approximation is rather poor for 
intermediate times. 

Because of the relationship between the MSD and the 
SISF given by Eq. [8j one may conclude that the stretched 
exponential form of I™(q,t) is an intrinsic property orig- 
inating from the polymer nature of the lipid molecules. 
For t values when G™ (r, t) is close to a Gaussian, the SISF 
for lipid atoms can be well approximated by Eq. ( pToj ) . 
This agree men t, combined with the behavior of the MSD 
from Sec. 



Ill 



and the small values of 72 (t) for all t, 
suggests that the decay of the SISF for the phospho- 
rus (P) atoms should be well described by a Kohlraush- 
William-Watt (KWW) function A(q) exp[-[t/r(q)]^]. 



FIG. 7: (Color online) SISF for the CT atoms for two scat- 
tering vectors qi — 1.42 A -1 and ^2 = 0.75 A -1 . The dashed 
curves represent the Gaussian approximation Eq. (10), the 
thin-solid curves the non-Gaussian approximation Eq. ([8]), 
and the thick (red) curves are the exact SISF calculated using 
Eq. 0. 



The parameters f3(q) and r(q) have a nontrivial q depen- 
dence. Larger q corresponds to smaller displacements 
and, hence, a shorter time scale. Smaller q probes dis- 
placements over larger length scales, thus it corresponds 
to dynamics on longer time scales. Fro m Eq. Jl0|) and 



f or the fits to ([Sr p (t)] 2 ) from Section III we expect that 



(3(q) « 0.6 for a range of q values, and we determine 
this range later in this section. For t > 30 ns ~ t\ the 
SISF crosses over to the simple exponential form char- 
acterized by /3(q) = 1. However, this Fickian diffusion 
regime can be observed only for small q values, which 
computationally are quite expensive to reach due to the 
long MD simulation time required to calculate the decay 
of the SISF. Thus, in most all-atom MD simulations the 
mainly explored dynamic region is the sub-diffusive one. 

To examine the q dependence of /3(q), and r(q) we fit 
the SISF calculated for the P atoms for t > 10 ps to 
A(q)exp[-[t/r(q)]PM]. Shown in Fig. § is /?( q) (circles 
and left axis). For q = 0.6 A -1 , we found (3(q) w 0.6 
(dashed line), which is in agreement with the exponent 
calculated from Eq. (i.e., the sub-diffusive exponent 
of the MSD). For smaller q values, (3(q) increases sharply 
towards the value (3 = 1. Note that the diffraction peak 
for the lipid head groups has been observed to be around 
q ~ 0.75 A -1 , which corresponds to the lipids being sep- 
arated by approximately 9 A. However, the stretched ex- 
ponential fits yield (3(q) = 1 only for q < 0.1 A -1 , which 
corresponds to a length scale of 62 A, around nine DMPC 
lipid diameters. Notice that (3(q) is approximately equal 
to the exponent for the sub-diffusive dynamics of the 
lipids for q values around the diffraction peak for the 
lipid head group. 

In the q — > limit, D(q) = l/[r(q)q 2 ] is equal to 
the lateral self-diffusion coefficient Dl obtained from the 
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FIG. 8: The q dependence of /3(q) (circles and left axis) and 
D(q) = l/[r(q)q 2 ] (squares and right axis) for the P atoms 
from fits of if \q,t) to exp[-[t/r(q)f {q) ]. The dashed line is 
f3 = 0.6 which was obtained from the MSD. The solid line 
is the diffusion coefficient Dz,-msd calculated from the linear 
region of the MSD for the phosphorus atom. 



asymptotic behavior of the MSD. D(q) is also shown in 
Fig. M (squares and right axis) from the results of the fits 
to I<~(q,t) for the P atoms. Shown as a solid line in Fig|8] 
is the value of the diffusion coefficient D^-msd given in 
Table [T] for the P atoms. It should be noted that, after 
a steep decrease with q, D(q) in Fig [8] seems to level off 
for q < 0.2 A -1 and is close to the value obtained from 
the linear fit to the MSD for the Phosphorus atom and 
the center of mass. 



C. Lipid Hydrogen Atoms 

According to Eq. ([6| the quasi-elastic scattering signal 
is modulated by the scattering cross-sections bi nc and 
b co h. To determine Dl one has to measure the inco- 
herent dynamic structure factor, which is dominated by 
the hydrogen atoms in the lipids. Therefore, we exam- 
ine lf{q,t) calculated for all the hydrogen atoms in the 
lipids. 

The light hydrogen (H) atoms follow the motion of the 
heavy atoms. However, being covalently bound to the 
heavy atoms in the lipid, the H-atoms also perform vi- 
brations and rotations about these atoms. Thus, one ex- 
pects that lj?(q,t) deviates (especially at shorter times) 
noticeably from I™(q,t) for the carbon, phosphorus, or 
oxygen atoms. Nevertheless, I^(q,t) should have sim- 
ilar characteristics in that stretched exponential relax- 
ation is expected for 0.4 < q < 2.5 and there should be 
a crossover to exponential relaxation for small enough q. 
Furthermore, we found that on ps time scales the dynam- 
ics of the lipid atoms depend on the location of the atom 
in the lipid . Thus, the motion of the hydrogen atoms 
will also depend on the location of the atom within the 



FIG. 9: The MSD calculated using all the hydrogen atoms in 
the lipids. The inset is a(t) calculated using Eq. |l]) and the 
dashed line is at the calculated value of /3 = 0.45. 



lipid in the ps time scale. In this section we examine 
what can be learned from the SISF about the dynamics 
of the hydrogens and the lipids. 

Shown in Fig. [9] is the MSD calculated using all the 
lipid hydrogens. The short time ballistic, the sub- 
diffusive region, and the diffusive region are clearly seen. 
Also shown in the figure is a H (t) given by Eq. 0, and 
we find (3 H = 0.45 (dashed line in the inset). It is impor- 
tant to realize that this value of /3 is much different than 
for the center of mass. 

is the direct calculation of lf(q,t) 
Gaussian approximation given by 
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Shown in Fig. 
(solid line), the 

Eq. ( flO] ) (dashed line), and the expansion given by Eq. ^ 
(red line) for wave- vectors around the diffraction peak for 
the lipid tails (q = 1.42 A -1 ) and for the lipid head group 
(q = 0.75 A -1 ). For the H-atoms the Gaussian approx- 
imation, Eq. ( [lQj ) ? is rather poor in the T2 ~ 0.01 ns to 
T\ rsj 10 ns range, but the first non-Gaussian correction, 
Eq. (|8|), is a very good approximation for all times. 

is I^(q, t) for several values of q. For 



Shown in Fig. 
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large q there appears to be two stages to the decay of 
I?(Qi namely an initial faster decay that is very nearly 
exponential, and a slower decay that is best described by 
a stretched exponential for t > 10 ps, i.e., when the sub- 
diffusive behavior begins in the MSD (see Fig. [9|. For 
this analysis we will focus on the slower decay. Again, 
we fit If(q,t), for t > 10 ps, to A(q) exp[-[t/r(q)]^] 
and the results for f3(q) and D(q) = l/[q 2 r(q)] are shown 
in Fig. [I2j 

For q = 1.4 A -1 (the peak value for the static structure 
factor for carbon atoms in the lipid tails) we found (3(q) « 
0.43, which is close to but lower than the value of f3 = 
0.45 obtained from the MSD through Eq. ([T]) (dashed 
line). For the distances corresponding to these q values, 
the lipid is undergoing sub-diffusive motion, Fig. |2j For 
q ~ 0.5, /3(q) rapidly increases to a value close to one, 
and D(q) appears to be leveling off, but our simulations 
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FIG. 10: (Color online) SISF for the lipid hydrogen H atoms 
for two scattering vectors q\ = 1.42A -1 and q2 — 0.75A -1 . 
The dashed curves represent the Gaussian approximation 
Eq. (10), the thin-solid curves the non-Gaussian approxima- 
tion Eq. ([8]), and the thick (red) curves are the exact SISF 
calculated using Eq. Q. 



FIG. 12: The stretching exponent /3(q) and D(q) from fits 
to the SISF calculated using the hydrogen atoms in the lipid. 
The dashed line is (3 — 0.45 which was determined from the 
MSD, see Fig. [5] The solid line is Dz,-msd for the center of 
mass. The inset is a comparison of D(q) for the phosphorus 
atoms (dashed line) and the hydrogen atoms (dashed dotted 
line) . 
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FIG. 11: The self intermediate scattering function calculated 
for the hydrogen atoms in the lipid for q — 2.1, 1.4, 1.0, 0.7, 
0.5, 0.3, 0.2, and 0.1 A" 1 listed from left to right. 

are to short to accurately obtain the asymptotic value. 
The solid line in the figure corresponds to Z)£ M -msd = 
1.62 x 10 -7 cm 2 s _1 obtained from a linear fit to the 
MSD for the center of mass. D(q) increases faster with 
q for the calculation using all the hydrogens than for the 
P atoms (see inset to Fig. [l2|, but, as expected, the two 
approach the same value for smaller q. For small D(q) 
is independent of the atoms used in the calculation, and 
it is from these values of q where the diffusion coefficient 
should be obtained. 

By looking at D(q) for the phosphorus and hydro- 
gen atoms a clear pattern emerges. The diffusion co- 
efficient obtained from the relaxation time would appear 
to be larger when determined from larger values of q and 



for shorter times. We examined the literature to deter- 
mine if this observation is compatible with experimental 
findings, and some indicative results are shown in Ta- 
ble |III| for l,2-Dipalmitoyl-sn-Glycero-3-Phosphocholine 
(DPPC). DMPC and DPPC are often used as model sys- 
tems. While DPPC has slightly larger acyl chains, as 
compared to DMPC (16 to 14 C atoms), our results for 
DMPC should be comparable to those for DPPC, show- 
ing the same trends. By comparing experiments that dif- 
fer by only the q range (experiments 2 and 3) we see that 
the larger Dl is found for the larger q. If we compare ex- 
periments that differ by basically the time scale involved 
(experiments 1 and 2) the associated Dl is larger for the 
shorter time scale. While the Dl listed are not compa- 
rable with our calculated D(q), the experiments are in 
qualitative agreement with our analysis. 



V. CONCLUSIONS 

In spite of extensive literature dealing with the experi- 
mental and computational study of the dynamics of lipid 
bilayers, the determination of the lateral self-diffusion co- 
efficient Dl of lipid molecules in the membrane remains a 
controversial problem. The results for Dl appear to de- 
pend significantly on the time and length scales probed 
by the employed measurement or computational method. 
This observation prompted researchers to propose dif- 
ferent models that have been predominantly influenced 
by models of diffusion in dense fluids. However, lipid 
molecules are polymers, characterized by flexibility and 
connectivity, and their lateral diffusion in the leaflets of 
lipid bilayers is qualitatively different from the diffusion 
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TABLE III: Diffusion coefficients for DPPC in the fluid phase determined by neutron scattering for different time and length 
scales. 



Experiment 


Technique 


T 


length scale 


time scale 


D 


i0 

2 m 

3 m 


QENS (IN5) 
QENS (IN10) 
QENS (IN10) 


45-55 °C 
41 °C 
41 °C 


0.2-1.6 A" 1 
0.07-2 A -1 
1-2 A- 1 


66 ps (63 fieV) 
5 ns (0.8 fieV) 
5 ns (0.8 fieV) 


3.5 x 10" 6 cm 2 /s 
1.8xl0" 7 cm 2 /s 
1.6xl0~ 6 cm 2 /s 



of molecules in dense fluids. This difference is reflected 
in the time dependence of the MSD. We observe a broad 
sub-diffusive region, not present in simple fluids, between 
the short time ballistic region and the long time Fickian 
diffusion region. 

Here we developed a phenomenological memory func- 
tion approach, which can be used to describe the behav- 
ior of the MSD of lipid atoms and molecules over the 
whole time interval spanning from the ballistic regime 
(t < 10 fs ~ rs) to the Fickian diffusion regime (t 
> 30 ns ~ ti). By fitting the mean square displace- 
ment, we were able to identify three clearly separated 
time scales that correspond to three different dynamic 
regimes in the lipid system. Overall, our memory func- 
tion fits to the MSD matched remarkably well ([5r(t)] 2 ) 
from the MD simulation. However, some features in the 
fairly broad crossover region between the ballistic and 
sub-diffusive regimes are missed by our approximation. 
Our memory function approach allows us to calculate 
the lateral self-diffusion coefficient Dl of lipid atoms and 
molecules through Eq. (J5|, and shows that, while the dy- 
namics of the different atoms are very different at ps time 
scales, the motion at t > 30 ns ~ n is best described by 
simple diffusion. 

We also investigated the consequence of the sub- 
diffusive behavior of the MSD by examining the first 
two terms in the cumulant expansion of the SISF. The 
first term is referred to as the Gaussian approxima- 
tion and is exact if the probability distribution G(r, t) 
of the displacements is Gaussian. By calculating the 
first non-Gaussian correction, characterized by the pa- 
rameter 72 (t) (Eq. [9J, we found that G(r, t) is nearly 
Gaussian outside the sub-diffusive region. While for 
the P atoms and the center of mass of individual lipid 
molecules G(r, t) remained nearly Gaussian even in the 
sub-diffusive region, G(r, t) for the carbon and hydrogen 
atoms showed significant deviation from Gaussian in this 
region. Interestingly, a Gaussian distribution for the cen- 
ter of mass, whose width did not increase as the square 
root of time, for lipid displacements has been observed 
in previous simulations [12] [14] . 

When G(r, t) is Gaussian, the SISF (e.g., measured 
in inelastic neutron scattering (INS) experiments of pro- 
tonated lipid membrane samples) can be expressed in 
terms of the MSD, Eq. (10). The sub-diffusive dynamics 



strongly influences the scattering from the lipids when- 
ever the probing scattering vector, and frequency, 
(j, correspond to lengths and times in the sub-diffusive 
regime of the lipid atoms and molecules. In this case 
the SISF can be fairly well fitted with a stretched expo- 
nential (KWW) function, A(q) exp[-[t/r(q)]^]. Note, 
however, that for a given q the stretched exponential will 
not be able to reproduce the short (ballistic) and long 
(Fickian diffusion) time regimes. We examined the q de- 
pendence of D(q) = l/q 2 r(q) and /3(q). For the phospho- 
rus atoms we find that (3(q) « 0.6 for q w 0.6 A , which 
agrees with the exponent determined from the MSD. 

We emphasize that in most INS experiments on lipid 
bilayers primarily one explores the sub-diffusive region 
and not the long time Fickian diffusion region, thus mak- 
ing the use of Eq. ( 12 ) inadequate. Therefore, it should 
come as no surprise that the values of the lipid lateral 
self-diffusion coefficient extracted from INS experiments 
using Eq. (12) have a strong q dependence and generally 



overestimate the real Dl obtained by other (e.g., FRAP) 
experiments. 

Furthermore, it would be interesting to investigate the 
temperature dependence of the analysis presented in this 
paper. In some liquids, it has been observed that ac- 
tivated events become more pronounced as the temper- 
ature is lowered. It would be reasonable to expect a 
similar behavior in lipid membranes, thus the diffusion 
could change character at lower temperatures and in the 
gel phase. It is also unknown how the picture presented 
in this work would change in the gel phase of the lipid 
membrane. Since the tails of the lipids are aligned, do 
the lipids behave as if they are more rigid? This may 
change the sub-diffusive region dramatically. The ad- 
dition of cholesterol or proteins to the lipids may also 
change the nature of the sub-diffusive region and dra- 
matically change the dynamics of the lipids. 
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